function ivp

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0        '
%  where y = (y1, y2, y3 , ..., yn) is an n-vector

% clear all previous variables and plots
clear *
clf

% time interval
tmax=20;

%  initial values
y10=1; y20=0; y30=1;  y40=0;
y0=[y10 y20 y30 y40];

%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],y0);  
%[t,y] = ode23s(@rhs,[0 tmax],y0); 

%
%  2-D plot:  plot solutions as function of t
%
hold on
plot(t,y(:,1),'r')
plot(t,y(:,3),'--b')

% commands to label each axes
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('y-axis','FontSize',14,'FontWeight','bold')

% command to put legend into plot
legend(' y_1',' y_2',1);

% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 

hold off

%  define f(t,y)
function dy=rhs(t,y)
dy=zeros(4,1);
ep=0.04; k12=1; c12=1;
dy(1)=y(2);
dy(2)=-y(1)-ep*( k12*(y(3)-y(1)) + c12*(y(4)-y(2)) );
dy(3)=y(4);
dy(4)=-y(3)-ep*(- k12*(y(3)-y(1)) - c12*(y(4)-y(2)) );













